The Anisotropy of MHD Alfvenic Turbulence 
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ABSTRACT 

We perform direct 3-dimensional numerical simulations for magnetohydrodynamic (MHD) 
turbulence in a periodic box of size 2ir threaded by strong uniform magnetic fields. We use a 
pseudo-spectral code with hyperviscosity and hyperdiffusivity to solve the incompressible MHD 
equations. We analyze the structure of the eddies as a function of scale. A straightforward 
calculation of anisotropy in wavevector space shows that the anisotropy is scale-independent. 
We discuss why this is not the true scaling law and how the curvature of large-scale magnetic 
fields affects the power spectrum and leads to the wrong conclusion. When we correct for this 
effect, we find that the anisotropy of eddies depends on their size: smaller eddies are more 
elongated than larger ones along local magnetic field lines. The results are consistent with the 
scaling law ku ~ k 2 ^ 3 proposed by Goldreich and Sridhar (1995, 1997). Here fen (and k±) are 
wavenumbers measured relative to the local magnetic field direction. However, we see some 
systematic deviations which may be a sign of limitations to the model, or our inability to fully 
resolve the inertial range of turbulence in our simulations. 



Subject headings: ISM:general-MHD-turbulence 

1. Introduction 

Many astrophysical plasmas, including the in- 
terstellar medium and the solar wind, often show 
magnetic fields whose energy density is greater 
than or equal to the local kinetic energy den- 
sity. In these plasmas the magnetic fields play 
a dominant dynamical role, mediated by magne- 
tohydrodynamic (MHD) waves. In the incom- 
pressible limit, there are only two types of lin- 
ear modes: shear Alfven waves and pseudo Alfven 
waves. While these two modes have different po- 
larization directions, they have the same disper- 
sion relation and propagate along the magnetic 
field lines at the Alfven speed. Therefore, the non- 
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linear interactions of wave packets moving along 
the magnetic field lines at the Alfven speed deter- 
mine the dynamics of incompressible magnetized 
plasmas with a strong background field. In this 
paper, we study the anisotropy of the MHD tur- 
bulence in this regime. We will refer to this tur- 
bulence as incompressible Alfvenic turbulence. 

Nonlinear processes and the corresponding en- 
ergy spectrum of incompressible Alfvenic tur- 
bulence are still among the most controver- 
sial problems in MHD. Since the pioneering 
works of Iroshnikov (1963) and Kraichnan (1965), 
the Iroshnikov-Kraichnan (IK) theory has been 
widely accepted as a model for incompressible, 
highly conducting MHD turbulence. The IK 
theory predicts EM(k) ~ Ex{k) ~ fc -3 / 2 from 
a Kolmogorov-like dimensional analysis. Here, 
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EM{k) and Ex{k) are the magnetic and kinetic 
energy spectra respectively. In this framework, 
two counter-traveling eddies (i.e. Alfven wave 
packets) interact and transfer energy to smaller 
spatial scales only when they collide, as they 
move in opposite directions along the magnetic 
field lines. Since the duration of such a collision is 
shorter than the conventional eddy turnover time 
by a factor of t v (l)/tA{l), this collisional process 
is inefficient and the spectral energy transfer time 
as a function of scale I (= t cas (l)) increases by the 
same factor compared to the eddy turnover time 
(l/vi) in ordinary hydrodynamic turbulence. Here 
t v (l) — Ijvi and <a(0 = I/Va are eddy turnover 
time and Alfven time respectively, Va = B/y/inp, 
and B is the rms magnetic field strength. When 
the external field is strong, as assumed in IK the- 
ory, this quantity is usually set to So, the strength 
of the uniform background field. If the spectral 
energy cascade rate 

e~^L~i^2 (1) 

is scale-independent and EM(k) w £jf(fc), then 
we obtain the IK energy spectra. 

The IK theory assumes an isotropic distribution 
of energy in k-space. However, many researchers 
have argued that anisotropy is an important char- 
acteristic in MHD turbulence (for example, She- 
balin et al 1983, Montgomery and Matthaeus 
1995). This anisotropy results from the resonant 
conditions for 3-wave interactions (or 4-wave inter- 
actions, when 3-wave interactions are null). The 
resonant conditions for the 3-wave interactions are 



h.2 = 



(5) 



ki+k 2 = k 3 , 



(2) 
(3) 



where k's are wavevectors and w's are wave fre- 
quencies. The first condition can be viewed as 
momentum conservation and the second as energy 
conservation. Alfven waves satisfy the dispersion 
relation: ui = VA|fc|||, where fc|| is the component 
of wavevector parallel to the background magnetic 
field. Since only opposite-traveling wave pack- 
ets interact, ki and k 2 must have opposite signs. 
Then from equations (2) and (3), either k\\ i or k^ 2 
must be equal to 0. That is, zero frequency modes 
are essential for energy transfer. If k^ 2 = 0, we 
have 



(Shcbalin et al 1983). Therefore, in the wavevec- 
tor space, 3-wave interactions make energy cas- 
cade in directions perpendicular to the mean mag- 
netic field. Since the energy cascade is strictly per- 
pendicular to the mean magnetic field, the actives 
modes in wavevector space have a slab-like geom- 
etry with a constant width. The implication is 
that the nonlinear cascade of energy works against 
isotropy in k space. Furthermore, it is important 
to note that equations (2) and (3) are true only 
when wave amplitudes are constant. In reality, 
nonlinear interactions provide a natural broaden- 
ing mechanism, following the uncertainty relation, 
At ■ Alu ~ 1. In particular, if a wave has a fre- 
quency less than or comparable to the nonlinear 
interaction rate, it is effectively a zero frequency 
mode. 

Goldrcich and Sridhar (1995, 1997) showed that 
in the strong incompressible shear Alfvenic tur- 
bulence regime (i.e. ~ u>), these arguments 
lead to a new scaling law with a scale-dependent 
anisotropy. In this model smaller eddies are more 
elongated. Their arguments are based on the as- 
sumption of a critically balanced cascade, U^Va ~ 
kj_v t} where k^_ and fcii are wave numbers per- 
pendicular and parallel to the external dc field 2 . 
The argument given above for 3-wave interactions 
makes it clear that k± will tend to increase un- 
til it becomes important in the plasma dynam- 
ics. The assumption of strong nonlinearity implies 
that wave packets lose their identity after they 
travel one wavelength along the field lines. Conse- 
quently the eddy turnover time ((fc_i_v;) -1 ) is actu- 
ally the same as Alfvenic time ((ku Va) -1 ). In this 
model, the cascade time, t cas (l) can be determined 
without ambiguity: t cas w (fc ± u z ) -1 w (k\\ Va) -1 . 
Since the cascade time is comparable to the period 
of the Alfven wave, the 3-wave resonant condition 
can be violated according to the uncertainty rela- 
tion Aui ■ At <~ Va/cm • t cas <~ 1. The quantity fcii is 
the width of the active region in wavevector space. 
Finally the assumption of a scale-independent cas- 



2 Later in this paper, wc will use k^ and fey, instead of k± 
and fey, to represent their scaling relation in a slightly dif- 
ferent, yet wc believe equivalent, viewpoint. Here, k± and 
fey are wavenumbers with respect to the direction of the 
local magnetic field, not the external field. 
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(4) 



2 



cade rate e ~ vf/t cas (l) ~ E waves V A /L gives 



fc|| <~ L 



2/3 r —1/3 / B waves 



Vl 



V A {k±L)-V* 



1/3 



1/3 



(6) 
(7) 



where E waves is the wave energy per mass. These 
formulae assume that all scales, from L on down, 
are within the inertial range of MHD turbulence. 
Here the typical k\\ should be interpreted as the 
size of the range of parallel wavevectors, corre- 
sponding to a given k±, that contain significant 
energy. 

Matthaeus et al. (1998) recently tested this 
model numerically, and showed that the anisotropy 
of low frequency MHD turbulence scales linearly 
with the ratio of perturbed and total magnetic 
field strength (b/B), a result which seems incon- 
sistent with Goldreich and Sridhar's model. To 
explain this scaling relation, they suggested that 
the region of Fourier space where the energy trans- 
fer takes place actively is given by 
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T~NL 



(8) 



where tnl is the eddy turnover time of the energy 
containing length L. Consequently, "the region 
of the wave number space where spectral transfer 
is most vigorous" has a slab-like geometry with a 
constant width proportional to 1/{tnlBq). 

All these theories (except the IK theory) share 
a common prediction for anisotropy: anisotropy 
should be more pronounced on smaller scales. 
Oughton et al. (1994) and Ghosh and Gold- 
stein (1997) already reported this scale-dependent 
anisotropy through numerical simulations. The 
former extended Shcbalin et al. (1983) 's 2-D cal- 
culations to 3-D cases. To measure anisotropy, 
they introduced the Shebalin angles, 6q, defined 
bytan 2 Q = (£ k\ |Q(k, t)| 2 )/(£ fcf|Q(k, t)| 2 ), 
where Q is vector potential A, magnetic field 
B, or current J, etc. Greater Oq means greater 
anisotropy. They found that 6a < Ob < 6j. If the 
energy spectrum of B scales as EM(k) cx k~ s , then 
the spectra of vector potential and current scale as 
EA(k) cx k~ s ~ 2 and Ej(k) cx fc~ s+2 , respectively. 
The spectra of vector potential has the steepest 
slope among the three spectra. This means that 
the vector potential is least strongly dependent 



on small scales (and the current is most strongly 
dependent on small scales). Therefore they con- 
cluded that anisotropics are more pronounced at 
smaller scales 3 . They also found a similar ordering 
for velocity field (and vorticity, etc). On the other 
hand, Ghosh and Goldstein (1997) calculated the 
Shcbalin angles MS c\ function of wavenumber bin. 
They found that the Hall-MHD 4 simulations show 
increased anisotropy at small scales (i.e. greater 
Shebalin angles at smaller scales). However their 
standard MHD simulations do not show increased 
anisotropics at small scales. We refrain from com- 
paring their work with ours because they used dif- 
ferent physics (the Hall effect) and their simula- 
tions are 2^ dimensional. The simulations given 
in this paper are 3-D standard MHD simulations. 



3 The width of the active region in Fourier space, fcy , is a 
function of k i . If 



i small k± 



large fcx 



then we will have the ordering of Shcbalin angles as ob- 
served by Shebalin et al. and Oughton ct al. However, if 
fey cx kj_, then we do not expect any ordering among the an- 
gles. In §3, we will find that, no matter what the true func- 
tional form feii = fey (kj_) is, Fourier transformation smooths 
out the true relation and leads to a fake linear relationship 
between fey and k±. This suggests that their results are 
in contradiction to our discussion in §3. However, the ap- 
parent linear relationship between fey and k± in §3 is not 
perfectly linear. Instead, we expect fey = cifcx + C2, where 
ci is a decreasing function of Bo and C2 depends on the fey 
of forcing terms (or initial values of fey for decaying turbu- 
lence). The presence of C2 does not seem to be important 
in our Fig. 3 in §3. However, it does affect the calcula- 
tion of Shebalin angles. That is, because of C2, the ratio 
k\\/k± = ci + C2/k± becomes a function of k±. At the 
largest energy containing eddy scale, kj_ ~ C2 and hence 
fey/fcx ~ ci + O(l) ~ O(l). But, at small scales, the ra- 
tio can be much smaller than unity. Therefore, because of 
C2, we can obtain a scale-dependent anisotropy: (fe||/fex) 
at small kj_ is greater than that at large fex- Note that 
this is a result of the initial conditions, rather than a true 
scaling relation. This will lead the ordering of the angles 
as observed by Shcbalin et al. and Oughton et al.. This 
interpretation is qualitatively consistant with their results. 
For example, the ratio fey/fex = ci + C2/k r ^ ak , therefore 
(tan#Q) -1 , approaches to a constant value C2/k p ^ ak as Bo 

becomes strong. Here k^ ak is the wavenumber of the peak 
of the energy spectrum. It is also possible to explain the 
increased anisotropy at high magnetic Reynolds numbers. 
If the magnetic Reynolds number increases, then k v ^ ak in- 
creases and, therefore, the ratio decreases. 
4 Hall MHD includes the Hall term, which is important at 
ion-cyclotron scales. In this paper, we consider only the 
standard MHD equations. 
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We note that none of the previous papers quan- 
titatively compared their results with particular 
theories of anisotropy. 

In this paper, we examine the scaling law for 
Alfvenic MHD turbulence numerically and resolve 
the controversy concerning the anisotropic struc- 
ture of the turbulence. Our results are consistent 
with Goldrcich and Sridhar's model. We describe 
our numerical methods in §2. In §3, we describe 
our results for anisotropy in wavevector space. In 
this section, we demonstrate that none of the scal- 
ing laws mentioned above agrees with our data, 
and explain why a straightforward evaluation of 
the distribution of spectral power does not cor- 
respond to a physically meaningful set of scaling 
laws. We attribute this to the systematic effects of 
large scale curvature of the magnetic fields. In §4, 
we determine the shape of individual eddies, avoid- 
ing the systematic error described above. We com- 
pare our results to Goldreich and Sridhar's model 
and give our conclusions in §5. 

2. Numerical Methods and General Re- 
sults 

We have employed a pseudospectral code to 
solve the incompressible MHD equations in a pe- 
riodic box of size 27r: 

— = (VxV)xV-Vj (VxB)xB+i/V 2 V+f+VP', 

(9) 

r)R 

— = B- VV- V- VB + r/V 2 B, (10) 

where f is random driving force and P' = P + V • 
V/2. Other variables have their usual meaning. 
Vao = B /(4:TTp) 1/2 is the Alfven velocity of the 
background field, which is set to be order unity 
in our simulations. In pseudo spectral methods, 
we calculate the temporal evolution of the equa- 
tions (9) and (10) in Fourier space. To obtain the 
Fourier components of nonlinear terms, we first 
calculate them in real space, and transform back 
into Fourier space. We use 21 forcing components 
with 2 < k < VT2. Each forcing component has 
correlation time of one. The peak of energy injec- 
tion occurs at k w 2.5. The amplitudes of the forc- 
ing components are tuned to ensure V s=s 1. Since 
we expect the perturbed magnetic field strength b 
to be comparable to or less than V (i.e. b < 1), 
the resulting turbulence can be regarded as strong 



MHD turbulence (i.e. V ~ Bq ~ B). The average 
hclicity in these simulations is not zero. However, 
our results are insensitive to the value of helicity, 
possibly because the box size is only slightly bigger 
than the energy injection scale. We use an appro- 
priate projection operator to calculate VP' term 
in Fourier space and also to enforce divergence- 
free condition (V • V = V • B = 0). We use 
up to 256 3 collocation points. We use integra- 
tion factor technique for kinetic and magnetic dis- 
sipation terms and leap-frog method for nonlinear 
terms. We eliminate 2Ai oscillation of the leap- 
frog method by using an appropriate average. At 
t = 0, magnetic field has only uniform component 
and velocity has a support between 2 < k < 4 in 
wavevector space. 

We use hyperviscosity and hyperdiffusivity for 
dissipation terms to maximize the inertial range. 
The only exception is Run 256P--B 1 ; where we 
use physical viscosity and diffusivity. The power 
of hyperviscosity is set to 8, so that the dissipation 
term in the above equation is replaced with 

-^ 8 (V 2 ) 8 V, (11) 

where we set the value of v% from the condition 
v h (N/2) 2h At w 0.5 (see Borue and Orszag 1996). 
Here At is the time step and N is the number 
of grids in each direction. We use exactly same 
expression for the magnetic dissipation term. We 
list parameters used for the simulations in Tabic 
1. We use the notation 256X-Y, where X = H 
or P refers to hyper- or physical viscosity; Y = 
BqI or BoO.5 refers to the strength of the external 
magnetic fields. 

Fig. 1 shows the evolution of V 2 and B 2 from 
Runs 256H-B 1 and REF2. Results of 256H-P 1 
are plotted as solid lines, and results of REF2 
are plotted as dotted lines. Both runs have sim- 
ilar results for the overlapped time, which means 
that the values of V 2 and B 2 are not sensitive 
to the spatial resolution. The magnetic energy 
density grows fast at the beginning of the simu- 
lations, as a result of the stretching of external 
field lines by turbulent motion. Since the external 
field is strong, magnetic forces soon become strong 
enough to balance the stretching effect. This bal- 
ance happens at t <~ 1 and, after a transient stage 
(1 < t < 5), the turbulence reaches the satura- 
tion stage (t > 5). In this case the saturation 
stage begins quite early, since the external uniform 
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magnetic field is strong and the magnetic diffusiv- 
ity is effectively 0. In general, when the exter- 
nal field is weaker and/or magnetic diffusivity is 
larger, the saturation stage occurs later. At the 
saturation stage, B 2 <~ 1.45 and V 2 <~ 0.6. Since 
b 2 = B 2 - B 2 (= 1) ?y 0.45, there is a rough energy 
equipartition between b and V. Note that, since 
V ~ B, the condition for strong MHD turbulence 
is met. 

3. Fourier Analysis: A False Scaling Law? 
3.1. Fourier Analysis 

In this subsection, we investigate the spectral 
structure of Alfvenic turbulence. It is important 
to note that we lose some information about in- 
dividual eddies when we perform a global trans- 
formation, such as the Fourier transformation. In 
particular, the scaling law given in this subsection 
may not be true for individual eddies (See §3.2 for 
details). 

We plot the 1-dimensional energy spectra 
Ex{k) and EM{k) in Fig. 2. Both spectra peak 
at the same k, which is also the scale of the energy 
injection. This reflects that there is a rough en- 
ergy equipartition between b and V at the largest 
energy containing scale. In general, Ex{k) always 
peaks at the energy injection scale. However, 
E]vi(k) peaks at a larger k than Ex{k) when the 
external field is weak. The inertial range exhibits 
of two different scaling ranges: a small k range 
where the slopes are steep and a large k range 
where the slopes are mild. We believe that the 
former reflects a true inertial range, whereas the 
latter is a result of the bottleneck effect. The 1/k 
bottleneck effect is a common feature in numeri- 
cal hydrodynamic simulations with hyperviscosity 
(see Borue and Orzag 1996, 1995). Interestingly 
enough, the slope of the bottleneck is actually 
steeper than —1. This might mean that the bot- 
tleneck effect is less serious in the MHD case. The 
kinetic and magnetic energy spectra have slightly 
different powers in the inertial range. Although 
the power indexes of both spectra are close to 
—5/3 in the (true) inertial range, the data are 
also compatible with IK theory, where the power 
index is —3/2. 

We plot normalized 3-dimensional energy spec- 
tra of Run 256H-i? l in Fig. 3. In the figure, we 



plot contours of same 

E 3 (k ± ,k ll )/E 3 (k ± ,0) (12) 

in the (k±,k\\) plane. As the energy cascades in 
the directions perpendicular to the mean field, this 
ratio drops as we move away from the k± axis. 
The figure implies most of the energy is confined 
in a region around the k± axis. The thickness of 
the region depends on the strength of the external 
magnetic fields (Fig. 4). We use the contours of 
Es(k±, k\\)/ Es(k±, 0) = 0.5 to measure the thick- 
ness. The angle 9o.s is the angle formed by the 
contours and the horizontal axis. We can see 
that tan 0o.5 (= fc||/fcj_) is proportional to b/B^. 
This result confirms the scaling relation found by 
Matthaeus et al. (1998). Both velocity and mag- 
netic fields have very similar structures. From the 
fact that fc|| < k± in the active region (see Fig. 3), 
one can conclude that eddies do have anisotropic 
structure: eddies are stretched along the direction 
of the mean field. Apparently, Fig. 3 suggests 
that ex k± and, hence, that anisotropy is scale- 
independent. No theory mentioned in §1 agrees 
with our result. However, it is not clear from the 
figure whether or not the true anisotropy is a func- 
tion of scale: although the contours show a linear 
relationship between k\\ and k±, the relationship 
doesn't mean that all individual eddies have the 
same major axis to minor axis ratio. As explained 
in next subsection, the rotation of eddies by large- 
scale waves in the magnetic fields can distort the 
actual scaling relation and lead to the linear rela- 
tionship shown in the figure. 

Fig. 5 shows tphase as a function of wavenum- 
ber k. This time scale is defined as the average 
correlation time At such that Fourier components 
at k = k have a phase shift of 60° with respect 
to the original phase. We plot the result for zero 
frequency modes (i.e. fc||=0 modes). We see that 
tphase ~ k" 1 (that is, fc^ 1 ). How can we interpret 
this result? If a turbulent structure, characterized 
by a wavenumber k, moves a distance I then the 
phase is shifted by roughly kl, even if the eddy is 
unaffected by the motion. This implies that a large 
scale velocity V will change the phase at a rate kV, 
so that tphase oc fc _1 . This implies that our calcu- 
lation of tphase is dominated by large scale motions 
rather than by the local cascade of energy, as long 
as the nonlinear cascade rate is proportional to k 
to some power less than one, which is generally the 
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case. This is an example of how large-scale fluctu- 
ations can complicate attempts to find physically 
meaningful scaling relations. 

3.2. Rotation Effect 

In the previous subsection, wc showed that 
anisotropy appears to be scale-independent. In 
this subsection, we will show that this apparent 
scaling relation is an artifact caused by the Fourier 
transformation. That is, we will demonstrate that 
large-scale modes in the magnetic field can distort 
the actual scaling law at smaller scales when we 
perform a Fourier transformation. In this man- 
ner, a straightforward evaluation of anisotropy 
in Fourier space is strongly contaminated by the 
curvature of the large-scale magnetic fields and 
does not reflect the actual local anisotropy when 
anisotropy is more pronounced at smaller scales. 
Consequently, figures 3 and 4 are compatible with 
any scaling law that predicts that smaller eddies 
are more elongated. 

Fig. 3 implies that eddies have anisotropic 
shapes: on average, eddies are stretched along 
the direction of B . However, not all eddies are 
aligned along the large scale field. The elongation 
of an eddy is determined by its interaction with 
the local magnetic field, not the background field. 
Since the large-scale magnetic field lines wander 
with respect to Bo, all smaller scale eddies have 
similar angular distributions around B . We illus- 
trate this effect in Fig. 6. 

In this way, we can explain the results of Gol- 
dreich and Sridhar (1995) and Matthaeus et al. 
(1998) simultaneously. Suppose that eddies are 
oriented according to the local field lines (Fig. 6a). 
Goldreich and Sridhar's result implies smaller ed- 
dies are relatively more elongated. When we per- 
form a Fourier transformation and measure the 
ratio of k\\/k±, what we actually measure is not 
the ratio of the minor axis to the major axis of in- 
dividual eddies. Instead, because the direction of 
the local magnetic field varies according to loca- 
tion and Fourier transformation is none other than 
a (weighted) averaging process, we actually mea- 
sure the ratio averaged over all possible orientation 
of the eddies (Fig. 6b) . The Fourier transforma- 
tion sees that L\ (and l\) is the minor axis and 
L 2 (and Z 2 ) is the major axis of an eddy. The ra- 
tio of Li/L 2 (= h/h) is determined by the degree 
of the wandering of the large scale magnetic field 



lines with respect to Bo and, therefore, the ratio 
is nearly constant for all eddies, regardless of their 
sizes and shapes. In fact, the ratio will depend on 
the tangent of the angle between Bo and B. Since 
cx 1/Li and fc|| cx l/£ 2 , the measured value of 
the ratio k\\/k± is nearly scale- independent. We 
expect the angle 9 (= 90° — 9 W ) between B and 
B to be 

tanO = b/B . (13) 
Therefore, we have 

sin6> = b/B = cos6 W7 (14) 

which is exactly the scaling relation found by 
Matthaeus et al. 

4. Measuring Eddy Shapes 

In this section, we analyze the shape of eddies 
in real space. As explained in the previous sec- 
tion, we assume that elongated eddies are aligned 
in the direction of the local magnetic fields. If we 
want to visualize the shape of eddies, we need to 
first identify the direction of the local magnetic 
fields. It is important to note that, although we 
use the term 'local magnetic fields' for simplicity, 
the fields are different from B(r). Let us consider 
an eddy of size I. The 'local magnetic field' of the 
eddy must act as the 'large-scale magnetic field' 
for the eddy. Therefore, the 'local magnetic fields' 
for eddies of size I , must be smoothly varying vec- 
tor fields whose characteristic length of variation 
is > I. Note also that another eddy of size I' (7^ I) 
at the same location can have a slightly different 
direction for the 'local magnetic field.' The 'local 
magnetic fields' are functions of position (r) and 
eddy size (I). Then, how can we define the direc- 
tion of the local magnetic fields? We implement 2 
independent numerical algorithms to calculate the 
direction of the local magnetic fields. 

In the first method, the local magnetic fields 
are calculated by 

B, = (B(n) + B(r 2 ))/2 (15) 

and the second order structure functions for v and 
b are given by 

F 2 "(i?,z)-<|V(r 1 )-V(r 2 )| 2 >, (16) 
F b 2 {R,z)=< |b(n)-b(r 2 )| 2 >, (17) 
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where R = |z x (r 2 — ri)|, z = z • (r 2 — ri) and 
z = B;/|B;|. That is, R and z are coordinates in a 
cylindrical coordinate system in which the z-axis 
is parallel to B; (Fig. 7). B(r) and b(r) are the 
total and perturbed magnetic fields at a point r. 
As usual, brackets denotes a spatial average. 

In Fig. 8 we plot the second order structure 
functions in z-R plane. The horizontal axis (z- 
axis) is parallel to B;. The contours reflect the 
shapes of the eddies. Eddies are clearly elongated 
along the local field lines, and smaller eddies are 
more elongated. The velocity and magnetic fields 
show different structure at large scales. However, 
their small scale structure is quite similar. 

In Fig. 9, we plot R-intercepts and z-intercepts 
of the contours. The R-intercept and z-intercept 
of a given contour can be regarded as a measure 
of k± and k» for the corresponding eddy scale 5 . 
Fitting the results for velocity fields gives 

f fcO- 69 , (256H - B o 0.5) 
w: All ~ I fc^' 70 , (256JT-B 1) (18) 
1 fco- 73 , (256P-B 1) 

On the other hand, for the magnetic fields we find 

f fcO- 64 , (256H - B o 0.5) 
b:k n ~l fco- 50 , (256JT-B 1) (19) 
{ fcO- 53 , (256P-B0I) 

The velocity fields show good agreement with the 
scaling relation, ku <~ k\ , proposed by Goldre- 
ich and Sridhar (1995). The power indices are 
insensitive to the strength of Bq or the form of 
viscosity (and diffusivity) . However, the magnetic 
field shows different scaling behavior. When the 
external field is moderately strong (256H-£?o0.5), 
the power index is very close to 2/3. On the other 
hand, for stronger external fields (256H-£?ol and 
256P-£>ol), the power indices are smaller than 2/3. 
The results are insensitive to the choice of viscosity 
(and diffusivity). It is not clear whether the ex- 
istence of a separate scaling law for the magnetic 
field represents a physical effect not included in 

5 We interpret fcy as the inverse of the major axis of eddies 
and k± as that of the minor axis. Since the major axis is 
assumed to be parallel to the local magnetic fields, fey is 
the parallel wavemimber with respect to the local magnetic 
field direction. On the other hand, in §3.1, ku is parallel 
to Bo. Therefore, fcy and k± in this section have different 
meanings from fcy and fc^ in §3.1. 



Goldreich and Sridhar's model, e.g. the first signs 
of small scale intermittency in the magnetic field 
distribution, or merely the failure of the numerical 
models used here to fully resolve the inertial range 
of strong MHD turbulence. 

In the second method, we employ a completely 
different approach. We obtain the local large scale 
magnetic fields by filtering out the small scale 
magnetic fields: 

B ff (r)=^B(r')0(|r-r'|), (20) 

where <j>(r) cx exp(— r 2 /af) is a gaussian function. 
To determine the shape of small scale eddies (i.e. 
eddies smaller than the filter size, ~ oy), we con- 
sider the following quantities: 

v CT (r) = V(r) - V CT (r), (21) 
b CT (r) = B(r) - B ff (r), (22) 

where and B CT are filtered fields (cf. eq. (20)). 
The fields v CT (r) and b CT (r) represent small scale 
fluctuation of velocity and magnetic fields, the 
shape of which can be regarded as an adequate 
approximation of small scale eddies. From these 
two fields we calculate the structure functions 

F v {R,z) = \v tr (r 2 )-v tr (r 1 )\, (23) 
F b (R,z) = |b CT (r 2 )-b ff (r 1 )|, (24) 

where R and z are similarly defined as in the first 
method with z || B CT (ri) (Fig. 10). 

In Fig. 11, we plot the results of the second 
method. We can clearly observe the flattening ef- 
fect: When the filter size is large, eddies are less 
anisotropic. When filter size is small, eddies show 
highly anisotropic structure. In the figure, we plot 
only the results for magnetic fields. Velocity fields 
show similar trends. The thick contours represent 
F b (R,z) =< |b CT | 2 >^ 2 . The values next to the 
thick contours are < |b CT | 2 > 1/>2 . 

This second method is useful for visualizing 
small scale eddies, but may not be as useful for 
quantitative analysis. The difficulty comes from 
the fact that there is no well defined eddy scale 
associated with filtered fields v CT (r) and b CT (r). 

In this paper, we will not pursue quantitative 
analysis using the method 2. Instead, we just wish 
to point out that both methods describe the same 
scaling law: smaller eddies are relatively more 
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stretched along the local magnetic field lines. In 
particular, the results from the first method are 
consistent with the relation k\\ ~ fc^ 3 proposed 
by Goldrcich and Sridhar (1995). 

5. Discussion and Conclusions 

Here we rederive the scaling law fey oc k 2 ^ 3 in 
the framework of 3-wave interactions. Except for 
the use of the uncertainty principle, the work in 
this section is independent of Goldreich and Srid- 
har's derivation. As noted by Goldreich and Srid- 
har (1995), 3-wave interactions are an adequate 
proxy for wave-wave interactions of all orders in 
a strong MHD turbulence. As long as we assume 
the locality of interactions, it is pointless to dis- 
tinguish k\\ and k± from k\\ and k±. Hence, for 
simplicity, we use k\\ and k± instead of k\\ and k± 
during the derivation. 

Suppose the 3-dimensional energy spectrum is 
given by 

£ 3 (fc±,fc||) = |V(k)| 2 cx^ 2 «.f, (25) 

where V(k) is the amplitude of the mode whose 
wavevector is k and f(u) is a positive, symmetric 
function of u (cf. equation (7) of Goldreich and 
Sridhar (1995)) which describes the power distri- 
bution as a function of eddy shape. If the width 
(or, thickness) of the energy spectrum in the di- 
rection of fc|| is kj_, then we can write 

£k(fc ± ,fc||)a*I 2a /(fc||/*£). (26) 

If the width is caused by the uncertainty principle 
(At ■ Aui w 1 with At oc t cas (l) and Alo oc fc||), 
then 

t cas (l) oc fc"' 3 . (27) 

Suppose the energy cascade rate e <~ vf/t cas (l), 
where I = 2n/k±, is scale- independent. Be- 
cause vf ~ k-^ a k\k^ (~ k±E(k±), E(k±) = 1- 
dimensional spectrum) and t cas (l) <~ fcj , we have 

fc- 2 «+ 2 + 2 ^ - const. (28) 

Therefore, 

l-a + = O. (29) 

Now, let's pick up a mode at a wavevector p and 
consider nonlinear interactions with other wave 



modes. First, the strength of the interaction with 
another mode at q is oc p|V(p)||V(q)|. (This 
comes from the (V x V) x V term in equation 
(9)) Hereafter we assume p = |p| « p_\_. Second, 
the number of interactions is oc p 2 p 13 . This is the 
number of modes near p. Here we use locality of 
3-wave interactions. If the interactions are ran- 
dom, the net change of amplitude per unit time 
will be the strength of the interaction times the 
square root of the number of interactions, or 

|AV(p)| ocp|V(p)||V(q)| • ( P V) 1/2 - (30) 

Therefore, we have 

t cas oc |V(p)|/|AV(p)| oc P -V P -W 2 , (31) 

where we assumed p oc q. Equating this with 
t cas oc p~^ ', we can write 

a -2 = -(3/2. (32) 

From equations (29) and (32), we have 

a = 5/3,/? = 2/3, (33) 

which is just the result of Goldreich and Sridhar 
(1995): 

fc,l oc fc 2/3 . (34) 

As a consequence, the 3-D energy spectrum be- 
comes 

£ 3 (fc±,fc||)oc^ 10/3 /(fc||A ± /3 ) (35) 
and the corresponding 1-D spectrum is given by 
E(k) oc k- 5 / 3 . (36) 

In summary, we have shown that the anisotropy 
of Alfvcnic turbulence depends on the spatial 
scales of eddies. In particular, our results confirm 
the claim by Goldreich and Sridhar (1995, 1997) 
that smaller eddies are relatively more elongated 
along the direction of the local magnetic field 
lines than larger ones. Quantitative measurements 
of the anisotropy using the velocity fields show 
good agreement with their proposed scaling law, 
fc|| ~ fcj_' as long we interpret these wavenum- 
bers as referring to the local magnetic field direc- 
tion. However, when the external magnetic field 
is very strong, magnetic fields scale somewhat dif- 
ferently, showing a slightly more rapid increase in 
anisotropy at smaller scales. 
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It is important to note that the correct scaling 
laws depend on comparing the eddy shape to the 
local magnetic field direction. 

As a final note, we wish to stress that our re- 
sults are not in agreement with the IK theory. The 
IK theory is based on the assumption of isotropy 
in wavenumber space, which may be true when the 
external magnetic field is very weak or zero. How- 
ever, even in these cases, the turbulence is globally 
isotropic, but locally very anisotropic. In this pa- 
per, we showed that eddies do show anisotropy 
and that the anisotropy is scale-dependent when 
there is a strong large scale field (which should 
apply to very small scales within any MHD tur- 
bulence cascade). On the other hand, our results 
are consistent with Goldreich and Sridhar's the- 
ory of strong MHD turbulence. More precisely, if 
we consider the ratio of hydrodynamic to Alfvenic 
rates, that is (feufe/fen Va), we find from equations 
(18) and (19) that 



k\V A 



Matthaeus, W.H., Oughton, S., Ghosh, S. and 
Hossain, M., 1998, Phy. Rev. Lett. 81, 2056 

Montgomery, D.C. and Matthaeus, W.H., 1995, 
ApJ, 447, 706 

Oughton, S., Priest, E.R. & Matthaeus, W.H. 
1994, J. Fluid Mcch., 280, 95 

Shcbalin J.V., Matthaeus, W.H. and Montgomery, 
D., 1983, J. Plasma Phys., 29, 525 



(37) 



where k w k±. From Fig. 2 we see that for the 
incrtial range this implies a ratio which is either 
constant or increasing with wavenumber. A con- 
stant ratio is predicted by Goldreich and Sridhar's 
model. The IK model predicts a slow decline. 

This work was partially supported by Na- 
tional Computational Science Alliance under 
CTS980010N and utilized the NCSA SGI/CRAY 
Origin2000. 
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Table 1 
Parameters 



Run a 


N 3 


V 




V 




B 


REF1 


144 3 


3.20 x 10" 


-28 


3.20 x 10" 


-28 


0.5 


REF2 


144 3 


3.20 x 10" 


-28 


3.20 x 10" 


-28 


1 


REF3 


144 3 


3.20 x 10" 


-28 


3.20 x 10" 


-28 


2 


REF4 


144 3 


3.20 x 10" 


-28 


3.20 x 10" 


-28 


3 


256H-£ 0.5 


256 3 


6.42 x 10" 


-32 


6.42 x 10" 


-32 


0.5 


256H-B 1 


256 3 


6.42 x 10" 


-32 


6.42 x 10" 


-32 


1 


256P-B 1 


256 3 


0.001 




0.001 




1 



a For 256 3 grids, we use the notation 256X-Y, where X = H or P refers to hyper- or physical viscosity; Y 
= B n l or £? 0.5 refers to the strength of the external magnetic fields. 
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Fig. 1 — Time evolution of V 2 and B 2 (Solid lines, 256H-B 1; Dotted lines, REF2). B Q = 1 for both runs. 

Fig. 2. — 256H-£> 1- Kinetic energy spectrum (E K (k)) and magnetic energy spectrum (E M (kj). For 
2 < k < 20, spectra are compatible with fc -5 / 3 . A 1/k bottleneck effect is observed before the dissipation 
cutoff k d ~ 90. 

Fig. 3. — 256H-_Bq1. 3D energy spectra normalized by the value on the k± axis. 

Fig. 4. — tan ©o,5 (= ku/k±) versus b/B n . These results are consistent with Matthaeus et al (1998) showing 
that anisotropy (fcii/fcjj scales linearly with b/B . 

Fig. 5. — 256H-£>ol- Phase correlation time t p hase as a function of wavenumber, defined as the average time 
required for the phase to shift by 60°. 

Fig. 6. Rotation effect, (a) Anisotropic wave packets are moving along magnetic field lines. We assume 
the wave packets are aligned with the direction of local magnetic field lines, (b) Large scale variations in 
the magnetic field direction causes wave packets to point in different directions as a function of location. 
We expect L\jLi = In Fourier space, L\ (and h) become <~ l/k± and L2 (and I2) becomes ~ l/fc||- 

Consequently, information about eddy shapes is lost when we look at the power spectrum. 

Fig. 7. — Method 1. We adopt a cylindrical coordinate system in which z-axis is parallel to B; = (B(ri) + 
B(r 2 ))/2. 

Fig. 8.- - 256H-£>ol- Visualization of eddies using method 1. Note that horizontal axis is z-axis (see Fig. 7 
for definition). Unit is grid spacing. Smaller eddies are more elongated. 

Fig. 9. — R-intercept (semi-minor axis; ~ l/fc±) versus z-intercept (semi-major axis; ~ V^ii) from Fig. 8. 

In 256H-_Bo0.5, both velocity and magnetic fields follow the relation fc|| <~ k^ . In 256H-_Bol and 256P-Sol, 
velocity fields follow the same scaling relation. However, magnetic fields scale slightly differently. 

Fig. 10. — Method 2. We adopt a cylindrical coordinate system in which the z-axis is parallel to large-scale 
fields B CT . 

Fig. 11. — 256H-£> 1- Visualization of eddies using method 2. Note that the horizontal axis is the z-axis 
(see Fig. 10 for definition). Unit is grid spacing. Only magnetic fields are shown. Smaller eddies are more 
flattened. 
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